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1.  Introduction 

The  authors  have  developed  a FORTRAN  program  that  can  be  used 
to  calculate  the  solution  to  the  homogeneous  system  of  first  order 
differential  equations  with  constant  coefficients  as  follows:* 


Y(t)  - AY (t) 

= Yo 


(1) 


Consider  Equation  (1)  and  let  Y(t)  be  an  n-vector  of  differentiable 
functions,  Yq  be  an  n-vector  of  constants  (initial  conditions),  and  A 

be  an  n X n matrix  with  constant  elements.  Then  the  solution  to 
Equation  (1)  can  be  written  as 


Y(t)  = $(t  - t0)YQ 


(2) 


where  $ is  an  n X n matrix  known  as  the  state  transition  matrix 
(fundamental  matrix)  whose  entries  are  functions  of  t.  For  example, 


i 

5 


| 

f 

j 


yx(t)  = y2(t) 
y2(t)  = -2y1(t)  - 2y2(t) 

is  a 2 X 2 system  of  linear  first  order  equations,  and  <I>(t) 
matrix 


(3) 

is  the 


e (cos  t - sin  t) 


e sin  t 


0 -t  . . 

■2  e sin  t 


-t 


e (cos  t - sin  t) 


(4) 


*The  symbol  A will  be  used  to  indicate  A is  a matrix  and  the  symbol  A 
will  be  used  to  indicate  A is  a column  vector. 
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To  determine  a solution  to  the  system,  given  a set  of  initial  conditions. 


one  need  only  premultiply  it  by  £(t  - tg) . This  result  has  been  in  the 

literature  for  quite  sometime^-  and  is  a special  case  of  the  more  general 
result  that 


Y(t)  = AY(t)  + BV(t) 

*<V  ■ ?0 


is  solved  by 


i(t  - t0)Y0  + 


t 

J !(t  - t0)BV(t)dt 

:o 


(6) 


(7) 


Clearly,  the  central  problem  in  determining  a particular  solution  in 

either  the  homogeneous  or  the  more  general  problem  is  the  calculation 

At  At 

of  $(t),  better  known  as  er-  . Subroutines  for  calculating  er-  have 

been  in  the  literature  for  quite  sometime.  All  those  known  to  the 

At 

authors,  however,  have  the  drawback  that  they  do  not  output  e—  in  a 

At 

form  similar  to  that  of  Equation  (4),  but  give  er-  a for  one  value  of 
t * tfi.  This  technique  is  widely  used  in  the  numerical  integration  of 

linear  systems.  Though  such  techniques  are  useful,  the  analytical  fo'.m 
of  the  solution  is  lost  together  with  time  constants  and  system  fre- 
quencies which  appear  in  an  analytical  representation  for  ^(t). 

The  computer  program  developed  in  this  report  cun  be  used  to 
obtain  £(t)  for  every  t.  The  user  need  only  input  the  system  matrix  A, 
the  dimension  of  A,  and  a set  of  error  tolerance  levels.  For  example, 
let 


*Frame,  J.  S.,  1 Matrix  Fyi.pc^9ps_and  Ajzpllc^tiOIVS,  IV >"  IEEE 
Spectrum.  June  1964,  pp.  123-131. 
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Then  the  program  output  will  contain  (among  other  things)  £(t)  expressed 
as  a sum  of  matrices  times  linearly  independent  functions: 


£(t)  = 


e C cos  t + 


-2  -1 


e t sin  t 


The  error  tolerance  levels  are  discussed  in  detail  in  Section  3. 

The  remainder  of  this  report  is  divided  into  three  sections. 
Section  2 gives  a description  of  the  problem  and  the  techniques  used  to 
solve  it.  Section  3 deals  exclusively  with  the  computer  program  and 
includes  a description  of  inputs  and  outputs.  Finally,  Section  4 gives 
a listing  of  the  program.  Those  familiar  with  analytic  functions  of 
matrices  can  go  directly  to  Section  3. 


2.  The  State  Transition  Matrix 

At 

a.  e-  . Definition,  and  Properties 
A t 

e—  , the  $(t)  of  Section  1,  is  defined  by  a power  series 
as 

* 2 2 3 3 

At  A t A t , _ 

e~  * I_  + At  + "21  + ”31  + • • • (10) 

A t 

where  £ represents  the  unit  matrix.  Since  e—  is  defined  analogously 
to  e01*",  a is  a scaler,  it  is  similar  in  many  respects.  For  example, 

A(t  + s)  At  As 
e—  * e—  er- 

d e^'t  1 At 
-JJ--A  e- 

e—  0 . I • (11) 


but,  generally  speaking 


(At  + Bt)  At  Bt 
e “ ~ t e”  c” 

At  Bt  Bt  At 
e~  e—  * e—  e— 


since  matrix  multiplication  is  not  a commutative  operation.  A compu- 

At 

tationally  useful  property  of  e—  is 


RAR  _ At  -1 

e Re—  R 


This  makes  it  possible  to  easily  calculate  e—  in  some  cases  by  trans- 
forming it  into  a similar  matrix  (e.g.,  normal  matrices  to  diagonal 
2 

matrices  ).  The  power  series  representation,  while  useful  for  numeri- 

At 

cal  work,  is  none  too  helpful  for  writing  e—  in  closed  form.  The 

At 

following  general  theorem  is  used  to  decompose  e—  into  a finite  sum 
of  n X n matrices  times  analytic  functions  in  one  variable. 

Theorem3 


If  f is  an  analytic  function  on  a simply  connected  open  set  D of 
the  complex  plane  that  contains  all  the  eigenvalues  A.  of  an  n X n 
matrix  B and  the  origin,  then  f(B.)  may  be  written  as  ^ 

“ "s  f <k-iv  > 

f®  * I I <k  - IK  *l.k  (14> 

j=l  k=l 


where  n^  denotes  the  multiplicity  of  the  j eigenvalue,  s is  the  number 

of  eigenvalues,  and  the  Z.  , are  n X n matrices  which  are  independent 

j 

of  f and  D (they  depend  only  on  the  matrix  1J)  . 

Since  e is  analytic  in  the  whole  complex  plane,  setting  B = At, 
f(x)  = eX  yields 

3 nJ  k-1  V 

At  V V tk  e i „ 


* -is 


(k  - 1)!  -j,k  • 


J-l  k-1 


Now  it  is  clear  that  e~  has  a representation  as  a finite  sum  of  n X n 
matrices  of  scalers  times  analytic  functions.  The  problem  is  to  deter- 
mine the  \.  and  the  Z.  . . 

j “j,k 


Berstein,  I.  N.,  Topics  in  Alaebra.  Blaisdell,  Waltham,  Massa- 
chusetts, 1964. 
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A review  of  the  terminology  involved  in  the  theorem  is  presented. 
A function  analytic  about  {0}  has  a power  series  representation  about 
{0}: 


f(z)  = ^ anzn  for  |z|  < 


therefore,  f(B)  can  be  formally  defined  by 


f(B>  - X • 


If  B is  an  n X n matrix,  an  eigenvalue  of  _B_  is  a scaler  \ . for  which 
a vector  x exists  such  that 


Bx  = \ .x 
J 


(B  - T)k  = 0 . 

J 


If  Equation  (19)  holds,  the  matrix  B - \I  is  singular,  and,  therefore, 

\ . is  a solution  of  the  equation, 

Det(B  - \I)  = 0 , _ (20) 

where  Det(B  - \I)  is  an  n*^  degree  polynomial  called  the  characteristic 
polynomial  of  B.  Its  roots,  which  are  just  the  eigenvalues  of  g.,  are 
also  called  the  characteristic  roots  of  B.  The  multiplicity  of  V is 
its  multiplicity  as  a root  of  Det(B  - \I) . ^ 

The  theorem  asserts  that  f(B)  exists  (the  power  series  of  Equation 
(17)  converges)  if  the  eigenvalues  are  in  D and  that  f(B)  has  a repre- 
sentation in  the  form  of  Equation  (14). 

b.  The  Characteristic  Polynomial  and  Its  Roots 

At 

Certainly,  the  first  problem  in  writing  e~  in  the  form  of 
Equation  (15)  is  to  calculate  the  eigenvalues  and  determine  their  multi- 
plicities. To  do  this  the  characteristic  polynomial  has  to  be  calcu- 
lated. If  , 

Det(A  - \I)  - \n  + d^\n  1 + . . . + dQ  (21) 


is  the  characteristic  polynomial,  the  coefficients  d^  are  calculated 
using  the  following  theorem. 


Theorem 


/ 


4 ■ tdki 


dA  - 1,  B.  - I»  B =0 
0 * -0  ‘ — n “ 


0 < k < n 


(22) 


Recall  that  tr(a^)  = £a^  *-s  c^e  sum  °f  c^e  diagonal  entries  of  (a^)  . 

Since  B^  * 0_  a good  way  to  check  for  errors  in  the  calculation  of  the 

coefficients  is  to  check  the  difference  between  the  calculated  B and 

—n 

the  theoretically  determined  values.  An  interesting  consequence  of 
Equation  (21)  is  that 

AB  . 

ZIlp  « I (23) 

n 


Ssi-A-i 

-dn 

if  A * exists  (d  * 0)  . 

— n 


(24) 


The  roots  of  Det(A  - \I)  can,  at  this  stage,  be  calculated  using 
any  one  of  a number  of  different  techniques.  In  this  program  the  classi- 
cal Newton -Raphs on  method  is  used.  The  roots  then  have  to  be  sorted 
and  multiplicities  counted.  Numerical  errors  can  be  generated  almost 
anywhere  in  our  program.  At  this  stage,  these  errors  necessitate  a 
decision.  For  example,  using  Newton-Raphson  the  equation 

x4  + 2x2  + 1 - 0 (25) 


will  have  calculated  roots  ± i,  e2  ± i where  e^,  e2  are  small 

(<  10  ^ using  a double  precision  version  of  an  IBM  root  routine),  but 
distinct  real  numbers.  The  solution  to  a 4 X 4 system  with  Equation  (25) 
as  its  characteristic  equation  will  be  calculated  to  be  of  the  form 


4 

Frame,  oj>.  cit. 
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z1>ie(6l+1)t  + Zj  .(er1),:  + Z3ile(€2+1)t  + Zlj4e(V1>t  . (26) 

But  since  the  roots  of  Equation  (25)  are  really  ±i,  each  of  multiplicity 
two,  the  solution  really  is 


— l,lelt  + — l,2telt  + — 2tle-it  + — 2,2te-it 


To  avoid  this  kind  of  problem  some  decision  has  to  be  made.  Namely,  a 
tolerance  e is  specified  such  that  a will  be  set  equal  to  b if 


a - b < e 


where  a and  b are  roots  of  the  characteristic  polynomial. 

For  example,  in  the  case  previously  considered,  e is  set  to  be 
large  enough  so  that 


[ej  - €2|  < e 


!•  & 


The  program  replaces  by  anc*>  rather  than  saying  that  the 

characteristic  polynomial  has  distinct  roots  ± i and  ± i>  claims 

that  it  has  a root  of  ± i of  multiplicity  two.  The  program  then 

checks  to  determine  if  the  root,  or  its  real  or  imaginary  part,  is  small 
enough  to  be  called  zero,  i.e.. 


a < e 


I ■ S 


If  Equation  (30)  is  satisifed,  a is  set  equal  to  zero.  In  the  example, 
the  final  output  is  two  roots,  ±i  (each  of  multiplicity  two),  if  e is 
sufficiently  large. 

Although  an  error  at  this  stage  will  make  the  solution  to  the  system 
look  radically  different,  it  will  probably  not  change  any  of  the  usual 
system  constants  in  a discontinuous  manner.  In  fact,  the  solution  to 
Equation  (1)  depends  continuously  on  the  eigenvalues.  To  see  this, 
notice  from  Equation  (10)  that  the  solution  depends  continuously  on  A. 
Using  Equation  (13)  it  may  be  assumed  that  A.  is  in  lower  triangular 

form 


Herstein,  loc.  cit. 
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(31) 


A * 


J?*  ent5*es  of  ^ are  just  the  eigenvalues.  If  the  roots  are 

^Six  Adhere  V th“  SyStem  iS  chanSed  to  one  with 


A'  = A + € 


(32) 


e = 


' e.  0 

\ 

\ 

t se 


Hence, 


A't  (A+e)t 
e—  = e ~ 


(33) 


which  is  continuous  in  the  e^. 

c • The  Constituent  Matrices 

The  matrices  are  called  the  constituent  matrices  and, 

as  was  previously  noted,  Ire  dependent  only  on  A,  not  on  the  analytic 
function  at  which  A is  evaluated.  Rather  than  launch  into  a general 

te^hn1ique  used  t0  calculate  the  constituent  matrices 
be  calcflated  for  a particular  example  and  the  illus- 
trated technique  generalized.  Suppose 


(34) 


with  characteristic  polynomial 

(x  1)  (x  ■ 2)  (35) 

which  has^roots ^ - 1 (of  multiplicity  two)  and  \2  = 2 (with  multiplic: 
one) . er-  must  be  of  the  form 


e%,l+te%.2+e2ti2>1 


(36 


10 
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Applying  Theorem  1 to  the  analytic  functions  f(x)  = x , g(x)  - x, 
2 

and  h(x)  = x and  substituting  A for  x,  the  resulting  equations  are 

1 = 1'-l,l  + 0 -1,2  + 1 -2,1 

— = 1— !,1  + 1—1,2  + 2— 2,1 


A - 1-Zlfl  + 2-Z12  + 4.Z21  , 

2 

By  using  matrix  notation  and  considering  I,  A,  A , £ ^ as  formal 
symbols  Equation  (37)  can  be  written  as  J * 


1 0 

1 1 2 )(  — 1 , 2 
1 2 m 


The  matrix 


10  1 
112 
1 2 4 


is  invertible  with  inverse 


*1.2 


or,  in  equation  form, 


0 2-1 


1 -2 


Z = 2A  - A 
i,i  ~ 


11 


(42) 


Zx  2 = -21+  3A  - A2 
Z2  i = I - 2A  + A2  . 


The  Z's  can  be  easily  calculated  from  Equation  (42).  In  general,  if  A 
is  an  n X n matrix,  the  column  vectors 


are  formed  and  yield  the  matrix  equation 


A 


a 


(45) 


where  V*”  is  the  transpose  of  an  n X n matrix  V,  called  the  Vandermonde 

of  Equation  (1).  The  simplest  way  to  calculate  the  entries  of  VJ  is  as 

t til 

follows.  Partition  V into  s,  n X n matrices.  The  k matrix  will 

be  filled  with  entries  calculated  from  the  k eigenvalue.  The  entries 
of  each  of  these  matrices  is  calculated  using  the  algorithm 


Vij-° 


ij  * 

V«-VI-lij-l+Vi.l#J 


i > j 
i * j 
i < j . 


(46) 


Ag  an  example,  suppose  A has  3 eigenvalues  of  multiplicity  3,  °f 
multiplicity  2,  and  with  multiplicity  1,  then  V1"  is 


12 


0 

0 

1 

0 

1 

1 

0 

X2 

l 

X3 

*1 

1 

A 

^2 

A 

3X1 

*1 

A 

3^2 

A 

4 A 2 

A 

*1 

io\3 

A 

5^2 

A 

t 6 

where  V is  always  invertible  so  Equation  (45)  can  always  be  solved, 
d.  Complex  Eigenvalues 

In  all  previous  examples  the  eigenvalues  have  been  real.  In 
general,  the  eigenvalues  may  be  complex  and  the  result  is  that  and 

Z may  in  general  be  complex.  The  program  developed  in  this  report 
i J 

handles  all  complex  computations  with  real  computation  and  no  FORTRAN 
complex  declarations  are  made.  As  a result,  double  precision  may  be 
used  to  provide  accurate  solutions. 

At 

If  the  eigenvalues  are  complex  then  the  matrix  e—  may  contain 
numbers  and  functions  which  are  complex  in  form  and  must  be  combined  to 
form  a solution  which  contains  only  real  numbers  and  real  functions. 

The  task  is  tedious  by  hand,  so  the  program  combines  complex  functions 
into  a real  form  convenient  to  the  user. 


The  generation  of  a real  form  for  e-  with  complex  eigenvalues  is 
handled  in  the  same  manner  as  in  the  case  of  a single  linear  equation 
with  constant  coefficients.  Namely,  if  \ = a + i£  is  a characteristic 

root  then  so  is  X.  = a - ig,  and  the  two  roots  have  the  same  multiplicity. 
The  term 


(k  - 1)! 


at  -igt 


-J,k 


h.J 


Frame,  loc.  cit. 
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At- 

occurs  in  . 


Write 


zjjk  = z+  izz 
“i,k  = — 

with  W,  WW,  Z,  ZZ  real  n x n matrices,  and  write  e±ipt 


(49) 

(50) 


as 


Multiplying  out  the 
becomes 


cos  0t  ± i Sin  0t 


(51) 


compxex  numbers  and  grouping  terms,  Equation  (48) 


..k-l  at 
t e 

(k  - 1)! 


cos0t(Z+ W)  + 


sin  0t(WW  - ZZ)  + i {sin  0t(Z  - W)  + cos  pt  (ZZ  + WW)} 


. At  . 

but  er~  is  reai  since  A is, 
must  vanish  for  all  t.  By 
connections 


(52) 


Lnereiore,  rne  complex  part  of  Equation  (c 
evaluating  Equation  (52)  at  0 and  rt/0,  the' 


ZZ  “ “HE  (53) 

are  established.  Simplifying  Equation  (52)  gives 
._k-l 

t at 

(k  - 1)!  e ^ cos  • Z.  " 2 sin  0t  • ZZ]  . (54) 

All  the  quantities  in  Equation  (54)  are  real. 

may  be  handled  using  ^nlj!  iMl^OT^UonsTScIir^at  comp^011  (45) 
rule?'3  “ ^ Ch“8ed  in‘°  2 * 2 «**«•  *-rdi„{ 


a + ip  /— 


r ■’) 

\P  a/ 


(55) 


If  z and  w are  complex  and  Z,  W are  their  corresponding 


matrices  then 


7,8 


Frame,  o£.  cit. 
^Herstein,  loc.  cit. 
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z + w 


Z + W 


Equations  (56)  are  sufficient  to  guarantee  that  algebraic  calculations 
done  with  the  matrices  will  agree  with  those  done  with  their  complex 
numbers.  This  monomorphism  is  extended  to  n X n matrices  in  the  natural 
way.  If  Z_  is  a complex  n X n matrix  and 


Z = A + IB  , 


(;  :) 


a 2n  X 2n  matrix. 


Similarily,  if  W,  C,  and  D are  n X 1 matrices  such  that 


W = C + iD  , (59) 

then 

"Ms)  • <60) 

9,10 

Equations  analogous  to  Equation  (56)  hold.  Under  these  transfor- 

mations, Equation  (59)  becomes 


Frame,  loc.  cit. 

Bradon,  G.  E.,  Introduction  to  Compact  Transformation  Groups. 
Academic  Press,  New  York,  New  York,  1972, 
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is  inverted  and  the  calculations  are  made  as  before. 

The  2n  X 2n  matrix  given  in  Equation  (62)  cannot  be  inverted  by 
inverting  the  two  n X n matrices  VR  and  VC  as  can  be  shown  by  this 
counter-example : 

(010 
0 0 1 
0-10 

hence 

Det(A  - \I)  = - (\3  + \)  (6' 

so 

\ - 0,  ±i  . (6 


The  real  matrix  corresponding  to  V*"  is 


n 


i 


i 


0 


1 1 

0 0 0 | 0 

0-1-1  I 0 0 


0 0 \ 

-1  1 \ 


0 


0 0 0 I 1 1 1 

0 1 -1  j 

0 0 0 J 


0 0 0 / 

0 -1  -1  / 


(67) 


None  of  the  3x3  matrices  are  invertible, 
e.  Error  Checks 


Recall  that 


At 


satisfies 


.At  . 
de—  . At 

dt  ' - 


(68) 


Furthermore,  among  all  the  analytic  functions  satisfying  Equation  (68) 
it  alone  satisfies 


f(0)  = I 


(69) 


Equations  (68)  and  (69)  can  be  used  to  construct  a test  on  the  validity 

At 

of  any  technique  purporting  to  calculate  er~  in  terms  of  its  constituent 
A t 

matrices.  If  e—  is  differentiated  with  respect  to  t [Equation  (15)] , 
Equation  (70)  results 


V V 

Z Z (k  - d:  [-j,kfi + xj-j,k] 


(70) 


j=l  k=l 

(without  loss  of  generality  we  can  set  Z.,n  +1  = 0).  In  terms  of 


Equation  (68) 


j’ 


s nj 


X t 


1 1 (T^irr  [hMi + ® ^7» 


j=l  k=l 


or  setting 


= Z! 


(k  - 1)!  -j.k  -j,k 


(72) 


n . 

s J 


X.t 


2 2tk'le3  + -a  • <73> 

j=l  k=l 

, . X.t 
lc“  1 i 

The  functions  t e are  linearly  independent  so  their  coefficient 
must  be  zero  for  the  left  hand  side  of  Equation  (73)  to  equal  the  null 
matrix: 


kZ ! . , + X .Z ! - AZ'  =0  k < n, 

“J.krfl  J j,k  — J jk  ~ j 


X.Z'.  - AZ.  =0  . 

~J,n.  - 

If  X.  = a + i0  is  complex  then,  using  analogous  notation. 


(74) 


kZR!  + cczr:  + gzi!  , - azr:  . = o 

— j.kfl  — J,k  j,k  j ,k  - 

kZI'.  ...  + azi'.  . - 3ZR!  . - AZI!  , = 0 

— J.kfl  — j,k  — j,k  j,k  - 


aZR!  + 0ZI'  - AZR!  - 0 

-J,n.  -j,n.  — J,nj  - 

aZl!  - 0ZR!  - AZI!  =0  . 

J »nj  J »nj  - 


(75) 


Checking  these  equations  provides  a good  overall  check  for  numerical 
errors.  The  program  keeps  track  of  the  maximum  entry  of  the  matrix  on  the  left 
hand  sides  of  Equations  (74)  and  (75)  (absolute  value  of  the  matrix  entry)  . 
After  all  calculations  have  been  done,  it  prints  out  this  maximum. 


3.  Program  Description 

a.  Flow  Chart  of  Program 

A flow  chart  of  the  program  containing  the  computations 
defined  in  Section  2 is  given  in  Figure  1.  All  blocks  in  the  flow  chart 
with  the  exception  of  the  last,  are  either  self-explanatory  or  have 
already  been  explained.  In  calculating  the  inverse  of  a matrix  the  IBM 

subroutine  INVERT, ^ which  employs  pivotal  condensation,  was  used.  In 


^Sy8tem/360  Scientific  Subroutine  Package.  International  Business 
Machines,  While  Plains,  New  York,  1968. 
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calculating  the  roots  of  the  characteristic  polynomial  POLRT,  another 

12 

IBM  routine  was  used.  It  is  a 500-step  Newton-Raphson  method  in  two 
variables . 

b.  Inputs  to  Program 

This  section  defines  the  form  of  the  input  data  for  the 
program.  The  firjt  four  cards  define  the  tolerance  constants  described 
in  Section  2.  These  cards  are  read  only  once  for  any  set  of  A matrices 
to  be  run.  If  the  dimension  card  (LDELi)  for  an  A matrix  is  zero  in 
value,  the  computer  run  will  terminate. 

The  order  and  format  of  the  input  data  cards  are  shown  in  Figure  2. 
The  variable  names  in  Figure  2 can  be  identified  as  follows: 

DB  - In  sorting  roots,  two  roots  which  differ  by  an  amount  less 
than  DB  are  set  equal.  This  tolerance  will  effect  the  multiplicity 
calculated  for  a root  and  hence  the  form  of  the  solution.  If  the  user 
selects  DB  = 0,  distinct  roots  are  likely  to  be  generated  changing  the 
form  of  the  solution.  Since  the  solution  to  Equation  (1)  depends 
continuously  on  the  eigenvalues,  the  form  of  the  solution  chosen  by 
the  program  will  give  the  correct  answer  in  a numeric  sense. 

DD  - If  the  determinant  of  the  Vandermonde  is  less  than  DD  in 
absolute  value,  an  error  message  is  printed  stating  that  the  Vandermonde 
matrix  is  singular.  When  this  message  is  encountered,  one  of  the 
following  occurred:  (1)  DD  was  selected  too  large,  (2)  User  entered 

data  improperly,  or  (3)  numerical  roundoff  or  truncation  error  is  the 
source  of  the  problem.  After  the  error  message  is  printed,  execution 
is  halted  and  the  next  data  set  is  read  in  (starting  with  LDEM) . If 
the  user  sets  DD  » 0,  the  program  halts  only  if  the  determinant  of  V 
is  zero. 

LDEM  - Order  of  the  system  matrix. 

A(I,J)  - Element  of  the  system  matrix. 


Sample  Program  Output 


Along  with  the  error  messages  already  discussed,  the  program 
will  print  the  A matrix,  its  characteristic  polynomial,  a table  of 

At 

characteristic  roots  and  their  multiplicities,  e~  expressed  as  a sum 
of  constituent  matrices  and  analytic  functions,  and  the  maximum  error 
as  determined  by  the  method  of  Section  2.e. 
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The  input  data  set  shown  in  Figure  3 was  used  to  generate  the 
output  shown  in  Figure  4. 


4.  Program  Listing 

The  listing  given  in  Figure  5 is  the  double  precision  version 
of  the  EAT  program. 


5.  Conclusions 

Programs  to  calculate  the  state  transition  matrix  (e^)  as  a 
linear  combination  of  functions  of  time  are  not  generally  available. 

At 

The  program  described  in  the  report  computes  e—  and  should  be  useful 
to  those  involved  in  the  solution  of  linear  differential  equations 
described  by  state  space  equations. 
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Figure  1.  Flow  chart  of  program. 
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CARD 

NO. 


VARIABLE  NAME 


FORMAT 


TOLERANCE 

1 

DB 

D15.0 

LEVELS 

1 2 

DC 

D15.0 

i 3 

DD 

D15.0 

FIRST  DATA 

4 

LDEM 

12 

SET 

< 

5 

A(1,1),  A(1,2),  A(1,3),  A(1,4) 

4D15.0 

SECOND  DATA 

i • 

' M 

LDEM 

12 

SET 

> 

M+1 

i 

A(1.1).  A(1,2),  AI1.3),  A(  1 ,4) 

4D15.B 

RUN 

1 . 

( K 

BLANK  CARD 

12 

TERMINATION 

\ 

\ 

Figure  2.  Input  data  card  format. 


( ? 

15 

1 

30 

i 

45 

i 

60  CARD  COLUMN  NUMBER 

0.0001 

0.0001 

0.0001 

0.0001 

4 

0.0 

1.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

1.0 

0.0 

0.0 

1.0 

1ST  DATA  SET 

-1.0 

0.0 

-2.0 

0.0 

3 

0.0 

1.0 

0.0 

0.0  1 

0.0 

1.0 

-4.0 

-6.0  | 

2ND  DATA  SET 

-4.0 

4 

0.0 

1.0 

0.0 

0.0  1 

0.0 

0.0 

00 

0.0 

1.0 

0.0 

0.0 

1.0 

3RD  DATA  SET 

-4.0 

-8.0 

-8.0 

-4.0 

4 

0.0 

0.0 

1.0 

00  ' 

1.0 

0.0 

0.0 

1.0 

4TH  DATA  SET 

0.0 

0.0 

0.0 

-1.0 

0.0 

1.0 

0.0 

0.0 

(BLANK  CARO) 

♦ F<  l)*Z<  It  1)  ♦ F<  Z) *ZZ ( 1,  i) 

♦ F<  3) *Z ( 1»  2)  ♦ F<  4)*ZZ(  It  Z) 


WHERE  f < I ) ARE 


Figure  4,  Sample  program  output. 


l 


— Ft  1)  = costx*  -.1OOOOOO0O*O1T 


Ft  3)  = < x**-  i )*COStX*  - . 1 00000000*011 


— FT “41  =“  TX*ir'l'1’*STMTX*~"-«10000000D'*OtT 


Awo  ntiERr  tht  z ma  zz  matrices  -are 


Z ( 1*1)  MATRIX -IS 


.1000000000*01  o.  - 0.  — — 0, 

0.  .1000000000*01  0.  0, 
0.  0.  .1000000000*01  0, 
0 . 0 . 0 . \ , 


ZZ(  1.  1)  MATRIX  is 

0.  . -.1500000000*01  0. 

.5000000000*00-  0.  -.5000000000*00  0 

0.  I .5000000000*00  0. 

.5000000000*011 -Oi .1500000000*01 0 


Z (1 . 2)  MATRIX  IS 

-.5000000000*00- 

.5000000000*00  0. 

O'.  .5000000001*00 

-.5000000000*00  0. 


.5000000000*00  0 

0 . - 

-.5000000000*00  0 


ZZ(  1.  2)  MATRIX  is 

-.5000000000*00  0. 

0.  -.51)00000000  *00 

.5000000000*00  0. 

— i ; -.5000000000*00 


MAXIMUM  ENTRY  OF  ERROR  MATRIX*  .407667520-12 


-.5000000000*00  0 

o. 

.5000000000*00  0 


♦ * * *-  * * #-#-  *-.«-*■*■*  • • 


1000000000*01 

500000000D*00 

5000000000*00 


5000000000*00 

5000000000*00 


5000000000*00 
.5000000000*00 


■ * * * * * * 


Figure  4.  (Continued) 


cac  Data  Set  Output 


i 

I 

| 


# 


A MATRIX  - 

0. 

■""0. 

-*<►000000000*01 


•1000000000*01 

0. 

-.6000000000*01 


0. 

• 1000000000*01 
-.<►000000000.01 


■*  *- 


* »-»  » 


* 'tt  • * 1 « 


CHARACTERISTIC  POLYNOMIAL 

( .^00000030*01)  ♦ 

< .600000000*01)  »x#*  1 ♦ 

I .^00000000*01 ) #x*»  2 ♦ 

( .100000000*01)  »X*»  3 


{ 

£ 

i 

j 


characteristic  roots 


Real  part  — - 
-.200000000*01 
-.looooooro+cr — 
-.100000000*01 


o. 


COMPLEX  part 

I 

.loooootroo+or 

.100000000*01 


MULTIPLICITY 

1 

1 

1 


♦ * * *-*  • 


EX»(AX)  = 

♦ E<  I ) *2 < 1,  1) 

♦ ?)*?(  2,  1)  ♦ P(  3) *22 ( 2*  1) 

WHERE  f<i>  are 

F<  11  * EXP(X»  -.200000000*01)  - 


Figure  4.  (Continued). 


i 

I 

4 

■I 


1 

] 
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F<  2)  = (X**  0)*EXP(X*  -.1000Q000D+Ql)*COS(X*  -.100000000+01) 
F(  3)  = (X**  0) *EXp(X*  -.1Q0000Q0D+01)*SIN(X»  -. 10000000D+01 ) 

ano  where  the  z ano  zz  matrices  are 


Z<  1«  1)  MATRIX  is 

•1000000000*01 

-.200000000D+01 

.4000000000+01 


.10000^000  + 01 
-.2000000000+01 
.4000000000+01 


.5000000000+00 

-.1000000000+01 

.2000000000+01 


Z<  2*  1)  MATRIX  is 

0. 

.2000000000+01 

-.4000000000+01 


-.1000000000+01 

.3000000000+01 

-.4000000000+01 


-.5000000000+00 

.1000000000+01 

-.1000000000+01 


ZZ(  2,  1)  MATRIA  is 


-.2000000000+01 

.2000000000+01 

0, 


-.2000000000+01 

.1000000000+01 

.2000000000+01 


-.5000000000+00 

0 . 

•1000000000+01 


MAXIMUM  ENTRY  OF  ERROR  MATRIX*  .605845180-27 


Figure  4,  (Continued) 


A *iAT4IX 

0. 

0.  

0. 

-.4000000000*01 


3rd.  Data.  Set  Output 


.1000000000*01 

o.  ■-  — 

0. 

-.8000000000*01 


.10000 00000*01 

0. 

-,8O0DDO0OOD*Ol 


0. 

0. 

.1000000000*01 

-.4000000000*01 


"CHARACTER  1 STIC  ‘TOL'r'MOHTAL 


f .400000000*017* 

■ r .'soooooooo^OTT_«nt»»'i  *~ 

( .8OOOOOOOO*Or)»7r**0  ♦ 

( .400000000*01)  «x*»  3 ♦ 

1 .roOOOOOO0*OtT_*X**'4  — 


CHARACTERISTIC  roots 


real  part 

-.100000000*01 

-.100000000*01 


COMPLEX  PART 
-.100000000*01 
.100000000*01 


MULTIPLICITY 

2 

2 


EXP (Ax ) * • 

♦ F<  l)*z<  1.  17  ♦ E ( 2)  *ZZ ( 1.1) 

♦ F(  3)*Z(  It  2)  ♦ F ( 4) *ZZ ( It  2) 


WHERE  FT!)  ARE 


Figure  4.  (Continued) . 
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f<  i)  = tx**  o>*EXPtx*  -.ioo©ooooo*on»cosrx*  -.io©ooo0ou*©i> 


"Ft  zr  * tx**  ot*E-xPTX*  «•.  iooooooou*©!  r*simx*  -riooooooop*otr- 


ft  j)  s (x»*  it*Exi»ix»  -.ioooooot>o*oi)»cosTx*-^.iut>t)t»(rooo*irn 


Ft  <►)  = tx**  it*exp { x * - . iooooooooto  1 ) * s i n < x* 


ano  where  the  z and  zz  matrices  xre 


Z(  1*  1)  MATRIX  is 


.1000000000*01 

0. 

0* 

0. 


.1000000000*01 


.1000000000*01 


or 

0. 

o r” - 

;iooooooooo*oi 


ZZ(  1,  1)  MATRIX  is 

-.2000000000*01 
.2000000000*01 
-.2000000000*01 
.4000000000*01  — 


30000000 0D*01 
.2000000000*01 
-.2000000000*01 
r 0000000000*01 


.1500000000*01 

.1000000000*01 

.2000000000*01 

.60oooooooo*or" 


-.5000000000*00 

.5000000000*00 

-.ioooOooooo*oi 
— ;20ooooooocr*ol 


Zt  If  21  MATRIX  IS 


-.1000000000*01 

.2000000000*01 

-.2000000000*01 


. 2000000 000*01 
.3000000000*01 
.2000000000*01 
.2000000000*01 


.2000000900*01 

,1000000000*01 

.2000000000*01 


.5000000000*00 

0.  

-.1000000000*01 


ZZ(  1,  2)  MATRIX  is 


-.1000000000*01 

0. 

.2000000000*01 

~-.4ooooooooo*ot 


-.1000000000*01 

-.1000000000*01 

.4000000000*01 

•V600000000D+01" 


-.5000000000*00 

*.iooooooooo*©i 

.3000000000*01 

■'.4000000000*01 


-.5000000009*00 

.1000000000*01 

"--.1000000000*01 


MAXIMUM  ENTRY  Of  ERROR  MATRIX®  . 100202Y30-10 


Figure  4.  (Continued). 


mm******** 


A MATRIX 

0. 


,IOOOOOOOOD*01 


0. 

0. 


4th  Data  Set  Output 


0. 

o.  - 

0. 

.1000000000*01 


.1000000000*01 

0. 

0. 

0. 


0. 

.1000000000*01 

-.1000000000*01 

0. 


characteristic  polynomial 

( .100000000*01)  ♦ 

<-o;  ) *x**  i ♦ 

( -.100000000*01)  *X**  2 ♦ 
(-0.  ) *X**  3 ♦ 

( .100000000*01)  *X*-*~4 


CHARACTERISTIC  ROOTS 


REAL  part 

complex  part 

multiplicity 

-.866035400*00 

.500000000*00 

1 

-.866035400*00 

-.500000000*00 

1 

,8660?5400*00 

-.500000000*00 

1 

.8660?5400*00 

.500000000*00 

1 

EXP (Ax ) = 


F! 

1)*Z< 

1* 

1) 

♦ F< 

2) *ZZ ( 1« 

1) 

F< 

3)»*< 

3* 

1) 

♦ F ( 

4)»ZZ<  3. 

1) 

Figure  4,  (Continued) . 
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tl 


MHERE£<I>  ARE 


Ft  1)  * (X«*  0)*EXP(X^  -.866O254O0*0O)*COS<X*  .SO0OOOOOD*OO> 
Ft  2)  a tX*«  0>*EXPtX*  -.86602540D*00)*SlNtX*  .50OOOOOOD*OO> 


— f-t-y  8"tl|K-»)  *£XPtX*-  .066025400*00 1 •COSfX'*  -.500000000+00) 


ft  « = tx**  ot*EXPTX*  .86602S4or>*oor*sroix*  -.500000000*00 j 


AND  KMERE  THE  Z AND  ZZ  MATRICES  ARE 


Zt  I#  1)  MATRIX  IS 


I 


.5000000000*00 
-.2A86751 350*00 
-.2886751350*00 

0. 


-.2886751350*00 

•500000000D*00 

0. 

-.2886751350*00 


-.5773502690*00 

0. 

.5000000000*00 

.2886751350*00 


2Z(  1.  1)  MATRIX  IS 


.2886751350*00 
' .5000000000*00 
-.5000000000*00 
-.5773502690*00 


-.500000000D*00 

^.2886751350*00 

.5773502690*00 

.5000000000*00 


-.5773502690*00 

.2836751350*00 

.5000000000*00 


rt  3.  n matrix  rs 


.5000000000*00 

.2886751350*00 

.2886751350*00 

0. 


.2886751350*00 

,500000000D*00 

0. 

.2886751350*00 


.5773502690*00 

0 . 

.5000000000*00 
-.2886751 35D*00 


ZZ < 3.  1)  MATRIX  is 


.2886751350*00 
-.5000000000*00 - 
,5000000000*00 
-.5773502690*00 


.5000000000*00 
-Y28 8 675 7350 *00 
.5773502690*00 
-.5000000000*00 


0 . 

- • 577350269  CffOtT 
.2886751350*00 
-.5000000000*00 


MAXIMUM  ENTRY  OF  ERROR  MATRIX*  . 379653230-28 


Figure  4.  (Concluded), 


o. 

-.5773502690*00 

.2886751350*00 

.5000000000*00 


,577350269D*00 

0. ” 

-.5000000000*00 
-. 2886751 35D* 00 


0. 

.577350269D*00 

-.2886751350*00 

.5000000000*00 


.5773502690*00 

o.  ■ 

.5000000000*00 

-.2886751350*00 
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PROGRAM  EAT  7A/74  OPT*I 


FTN  A. 2*75067  0S/16/7S  09.59.Z9. 


~ PROGRAM  EATI INPUTVOJTPUT  ,TAPE5*INPGT7TAPE6=0UTPUTJ 
0 1 HENS I ON  XCOF (10) »C0F (10) .ROOTR(IO) • ROOT I (XO). 

'DIMENSION  STORXTTOV.STORY(IO)  T~  • 

OIMENSION  VRUO.IO)  ,VC(I0.10)  ,RSAV(I0,I0) 

DIMENSION  KFUOr.NaOr.DELVUOl  .CMPXXUO) 

OIMENSION  A(IO.IO) ,Z<10.10) ,ZZ(10,I0) ,PRO0(20,2Q) 

OIMENSION  LI (20) .Ml (20) 

DOUBLE  PRECISION  aELL.CHPXX.A.Z.ZZTPROO 
DOUBLE  PRECISION  V,S,0,B 

DOUBLE  PRECISION  OA.W.TJCtDO  ' 

double  precision  xcof.cof.rootr.rooti.story.storx 
~ TJOUBLE" PRE  CTS  rON“Y«X 
OOUBLE  PRECISION  YY 
DOUBLE  PRECISION  YRtVCYRSAV 

READ  IN  ALL  THE  DATA" " 

REAO(S.lOl)  OB  

READ <S. 101)  DC 
READ(5.I01)  00 
ioi  format < 1015.0 > 

866  WRITE (6.7SSI 
WRITE (6,937) 

937  FORMAT!  601 1 XTTH»11 ' ~ 

WRITE (6,80 1 ) 

8oi  format  am  > 

READ (5, I)  LDEM 

K=0  - " 

IFF«0 

IFCLOEM)  75r.-T6r.-76I — 

761  REA0(5,88)  ((AU.J)  • J*I *LDEM)  , 1*1 , LDEM) 

88  F0RMAT(A0I5.0)  _ 

WHITE(b,60) 

60  FDRMAT(//////3XiBHA  MATRIX) 

CALL  WRITEULOEM.A) 

WRITE  (6,'7S5) ‘ 

WRITE<6,937) 

C 

C CALCULATE  THE  CHARACTERISTIC  POLYNOMIAL 

C 

XCOF ( LDEM. 1)« I. DO 

CALL"  TRACE  ( A.T0EM5DJ T 

XC0F(L0EMI»-8 

00  AO  I«I.LDEM  “ " ' 

00  AO  JM.LDEM 
AO  VR(I, J1«A(I, J) 

00  10  i-i.loem 

1 0 VR ( I , r> * VRTr,lT*TOJF  (LDEM)  

MM*L0EM-1 

M«LOEM  ' — ~ 

00  33  1*1 ,MM 

CALL  MATPRO(M,A,YR,PROO)  — 

call  tracE2<»roo«m,8) 

xcoF<M-i)*-ri.D0/o8Letrto«Trr»iiTT»9 

00  20  IP*I,M 

00  20  JP*1,N  
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— . - 20  vfniPi  JP)  «PROtnTP'."jpr 

DO  30  IP*I,'M  ' 

(SO  ' 30  VR ( IP. IP) *VR ( I P. IP)  ^CORM-TT  ’ ; ; 

33  -CONTINUE 

WRITE I6*75S)  

WRITE(6,104) 

. ... . F0RHATOX,?5HCH*HACTERrSTIC  P0LYN0MTAE77r 

65  CALL  WRITEP(XCOF.LOEM)  . _ . 

C CALCULATE  THE  CHARACTERISTIC  ROOTS 

CALL  POLRT(XCOF,COF,LDEM,ROOTR«ROOTI.IER> 

70  • 00  399  -TOVLOEW 

ir(ROOTR(I).LT.3d.ANO.ROOTR(I).GT.-Da>  ROOTR (I) *0.00 
IF(ROOTI(l).LT.aa.ANO.ROOTl(n.GT.4DBI  W>0TTaT*ff.T)0 

399  CONTINUE  ( 

IF  HER)  61,70*61  ' r 

7S  61  WRITE(6,I08>  1ER 

108  FORMAT  (///3X,"3DNERA0iT'IN'  ROOT  CALCIJL*TTO^HOOE',I2T 

GO  TO  866  _ 

C SORT  ANO  CLASSIFY  THE  ROOTS 

flO  C .............  ....... 

70  CALL  SORT (ROOTR, ROOTI.RELL. CMPXX. STORX. STORY, N.L0EM,K.D8) 

WRITE (6,755)  " 

75S  FORMAT (3X///) 

. ....  . WRITE(6,5S|  — ‘ 

S5  S5  FORMAT (2X.20HCHARACTERISTIC  ROOTS  //2X.9HREAL  PART.I3X.12HC0NPLEX 

1 PART,13X,12HHULiTlPLlClTY  ) 

IQ*K 

— — oo  roo  1*1,10 — ; 

700  WRIT£(6,701)  RELL(I>  .CMPXX (I ) ,N(I) 

oO  701  FORMAT  (1X«DIS,8,5X»015,B,I5X, 15)  * 

WRITE (6.755) 

WRITE (6,937)  - --  - - • - • - — =~ 

C 

— C — CALCULATE  1P«r  INVERT  THE  VANOCRHONPE 

95  C 

• KF(I)*I  - * - 

00  883  I«I,LOEM 

KFX«I*I  ’ ' - ' ' J 

B83  KF (I*I)*KFX*KF (I) 

•300  DO  313  III*I,K - . — 

IF (CMPXX ( 112 1 .NE.O.OO)  GO  TO  129 

313  CONTINUE  

CALL  VAN0R(R£LL,L0EM,N,K,V! 

GO  TO  107 

lf)5  129  L0EM2*2*L0EH 

CALL  VAN0(REtl»CNPXX»LOEN»N»K» V)  ' — * *' 

CALL  ARRAY (2,LOEH2,LOEM2,20, 20, S.V* 

CALL  OlNVRT(S,L0EH2,0,LI,MI)  " 

CALL  ARRAY ( I,LOEM2,LOEM2,20,20 ,S,V) 

110  U«0A8S(0>  

0«0S0RT(0) 

IFtO.OT.OO.OR.UiCT.^OO)  00  TO  5*6  ' " ~ * 

647  WRITE (6,646) 

648  FORMAT (3X.40HW  MATRIX  IS  SINGULAR  — CHECK  INPUT  0ATA1  
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PS  ' ' GO  TO  866  — — ; * 

696  KFG=1 

GO  TO  555  ■-  — 

107  CALL  ARRAY <2,LOEM,LOEM, 20*20. S.V) 

CALL  G1NVRT(5»L0EM,0,LI.M1)  ‘ 

120  CALL  ARRAY (1 tLDE4.LDEM.20  *20  ,S,9) 

lF(O.LT.T)O.AMO.O^GTr*OD)  “GO  “TOnbAT 

C 

C OETERMTNE  THE  FORM  OF  EXP(AT)  " 

C 

125  649  KFG*0  . - - 

LOEM2«LOEM 

555  1K*0 ■ “ " 

MRITE<6,756> 

756  FORMAT (3X///1  ~ - 

170  MRITE(6.4ll) 

IIK«0  . . 

KOOP«0 

00  444  I«I  :K  - — 

IF(KOOP)  414.415.414 

135  414  KOOP«0  - . . 

GO  TO  444 
415  LGM*N( I ) 

00  400  J«I,LGM 

IIK«IIK*1  — - • 

1 40  IF<CMPXX<I|>  406.401  ..06 

401  MRITE(6,410)  IIA.I.J 

GO  TO  400 

406  LOG*I IK*1  

WRI TE ( 6.440 ) IIK.I.J.LOG.I.J 

1 45  I IK*LCG  ' ' • “ •-  

KOOP*l 

400  CONTINUE 
444  CONTINUE 

411  format <3x.9he*P(AA)  » /i 

ISO  440  FORMAT <3X.4H*  F ( . I2.4H) *Z ( * 12* lH» . 1 2.6H)  ♦ F ( .I2.5M) *Z2 ( .12, 1H.. 

1I2.1HI/I  --  

410  FORMAf (3X.4H*  F ( . I2.4H) *2(>I2.1M,,12.1H)/) 

MRITE(6,4301 

KOOP*0 

155  IIMO 

DO  901  1*1,* 

LGM«N<I)  — 

IF(KOOP)  614.615,614 

614  KOOP*0 

1*,0  GO  TO  901 

615  00  600  J«I,LGM 
IIK-IIK.l 

IF(CMPXX(Ill  606,601,606  

601  IF(J-l)  604,602,604 
1 65  602  MRITF (6,603)  IH.RELLII) 

GO  TO  600 
604  JJJ»J-1 

MRITE(6,605>  1 1 A , JJJ.RELL 1 1) 

GO  TO  600  ' — - 

170  606  KOOP*l 

IF (RELL (III  612,607,612 
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. 60T-  fTt>tl-*l»*6*«*«0  - 

608  WRITE<6,609>  IK.CMPXX(l) 

I IK*l IK*1  - 

WRIT£(6,709)  1 1X.CMPXX < I ) ' 

KOOP=l  - _ - - — 

GO  TO  600 

....  610  JJJ3J-1  - --  --  • ~ 

WRITE (6,61 1 ) 1 IK , JJJ.CMPXX ( 2 ) 

I1K*11K*1 

WRITE (6.71 1 ) llK.JJJ.CHPXX(l) 

GO  TO  600  ‘ ~ 

612  JJJU-1 

WR1TET6.613)  IIK.  JJJ.RELLd)  .CMPXX  ( l > ' " ' 

UK*I1K*1 

WRITE <6. 713)  IIK.JJJ.RELL ( 1 > .CMPXX ( I > 

600  CONTINUE 
901  CONTINUE 

605  FORMAT (///3X.2HFI .12.8H)  - (X»».12,9H)»EXP(X»»015.8.1H) > 
6lT  FORMA TT///3X.2HF( ,12.870  * <X»», 12. 8H) »COS (X».0l5.8. IH>T 
711  FORMAT (///3X .2HF ( , 12.8H)  * (X»».i2,9M)*SlN(X»»015.8.1H) ) 

613  FORMAT(///3X.2HF( .12.8H)  = (X»« .12, 8H) »EXP (X». 015.8. 
18H)»COS(X», 015.8. 1H)) 

713  FORMAT ( ///3X.2HF ( . I2.8H)  « (X»». 12. 8H) «EXP (X*. 015.8, 
18h)«S1N(X«, 015.8. 1H» ) 

603  FORMAT(///3X.2HF(,I2.10H)  * EXP(X», 015.8, 1H))  

609  FORM*T(///3X.2HF( .12.4H)  « .6HC0S ( X».Ol 5.8, 1H> ) 

709  FORMAT  (///3X.2HFI,  1 2I54H)  * ,6HS1N(X».01S.8,1H) ) 

430  FORMAT (///3X, l 4HWMERE  F<1>  ARE) 

WRITE(6,760i 

760  FORMAT (///3X.35HANO  WHERE  THE  Z ANn  n MATRICES  are  >_ 

C CALCULATE  THE  CONSTITUENT  MATRICES 
C 

K00P«U 

lK-0 

YY»0.U0 

—00  809  1«1,X  * - 

iF(KOOP)  616.617,616 

616  KOOP«0 
1K«1K*1 
GO  TO  809 

617  MMM«N(I) 

00  109  J»1,MMM 
tK«tK*l 

IF(CMPXX(1))618,619,6I8 

618  KOOP«l 
KFG«1 

GO  TO  620 
— 619  KTG*0 

620  00  111  L1L*1 ,LJEM 
00  til  LJL«I.LOEM 
lF(LIL-LJL)  113,112.113 

112  Z(L1L,L.IL)«V(1K,1) 

GO  TO  111 

113  2n.tL»L.t)«0.0 
111  CONTINUE 

!F(KFO>  205,206,205 


Figure  5.  (Continued). 
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-205-00  11#  lttslTtOC* ~ 

00  114  JLL*1.L0EN 

IF  < ILL* JLL ) 116.115.116  

115  MKMsIK*LO£M 

ZZ<ItL*JLLl«V(MKM.11 

GO  TO  114 

H*  rzrittrJttr  *«vo -■ 

114  CONTINUE 

206  «0  t48'-Jt*2.t0EN  " - - 

IF(JL.EQ.2>  SO  TO  173 

• IFTJt.E0.3T  SO  TO  132 

CALL  MATPROILOEN.A.RSAV.PROO) 

ccrT0  273-  - — 

132  00  236  lPM.lOEM 
DD  236  JP»1.LDEN 
236  RSAV(1P.JP)»A(1P.JP> 

CALL  MATPRO(LDEM.A.RSAV.PROO) 

273  00  199  lPM.LOe* 

00  199  JPM.LOEN 

199  RSAVIIP. JP) »PROO< IP. JP) 

GO  TO  120 

173  00  -919  1P>1«L0EM 
00  919  JP»1.L0CN 
919  PR00(IP.JPl«A(IP.JP) 

120  00  119  tN«l,LOt*  ' ' 

DO  119  JNM.LKN 

Z(1N.JN)*Z(IN«JN).V(1K*JL) •PRDOTIN. JN> ‘ 

IF(KFG)  333.119.313 
333  LKL«1K*L0CM 

ZZ(1N.JN)»ZZ(1N,JN).V(L3L.JU*PH00(IN.JN) 

119  CONTINUE - 

140  CONTINUE 

IF(KF'i)  758.757.758 
758  00  759  N1N-1.L0E8 

00  759  NJN.l.LQCN  ■ -- 

Z(N1N.NJN»«2.00»Z(NIN,NJN) 

• ~ 759  ZZ(NlR.NJN)»2.Tra«ZZ(NTN,NJN)' 

757  8PITE (6.778)  I.J 
IF(J.EO.l)  X«l.DO 
if (J.gt.1)  a*08le (Float (aF(j-i ) ) j 
00  827  4U«1.lO£4 
00  827  MUT.l.LOEM 

2OW,WTl«7(MUi50T)/X — — ' 

827  CONTINUE 

00  321  NU»1.L0CM 
00  121  NUT.l.LOE* 

*«Z(NU.NUT) 

IF(X.LT.OC.ANO.X.ST.-OC)  Z(MU.MUT)>0.00 

321  CONTINUE  - • - 

CALL  KN1TEI  (L0E9.Z) 

IF(KFG)  207.309.287 
109  AaRELLd ) 

IF ( J.EO.l ) 50  TO  15T 
CALL  NATPRO ( LOCK. A. 9R. PROD! 

do-os?  ip-t,Loe* — — 

00  037  JP^I.LOEN 

T«X*»R(IP.  JP)  *Z(1P.  JP)  »f  LOAT  ( J-l ) -»*00(  IP.JPt 


Figure  5.  (Continued) , 
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Y»0A8S(Y) 

14  IF(Y-YY)  837.17.17 
17  YY*Y 

837  CONTINUE 
157  00  783  IP=1.L0EN 
00  783  JP«1,L0EM 
783  VR(IP.JP)«Z<IP.JP) 

IF(U.NE.NHM)  GO  70  109 
CALL  NATPRO(LOEN.A.VR.PROO) 

00  387  IP»l,L0EN 
00  387  JP»1.L0EN 
Y»PROO(IP.JP)-¥RUP.JP)»X 
Y*0A8S(Y) 

15  IF(Y-YY)  387.16.16 

16  YY*Y 
387  CONTINUE 

8S6  IF(KFG)  207.109.207 
207  8RITE (6.779)  I.J 
IF(J.EO.I)  X-I.OO 
IF (J.GT.l)  X*0BLE  (FLOAT  <XFU-1))  > 

00  839  HUM.LOEN 
UO  839  MUT.l.LOEN 
22 (NU.NUT ) *ZZ  (HU.NUT ) /X 
839  CONTINUE 

00  128  NUM.LOEN 
UO  128  HUT.l.LDEH 
XaZZINU. NUT) 

IF(X.LT.OC.ANO.X.GT.-OC)  ZZ (HU.NUT) *0.00 
128  CONTINUE 

CALL  2RITE1 (LOEN.ZZ ) 

IF (J. EO. I)  GO  TO  717 
CALL  NATPRO (LOEN.A.9R.PR30) 

00  703  IP*1,lOEN 
UO  703  JP*1.lOEN 

Y«FLOAT(J-I)»Z(I?,JP).R£LL<I»»¥R(IP. JP> *CNPXX (I ) »¥C (IP. JP) 
Y*Y-PR00(IP.JP) 

Y«OABS(Y) 

21  IF(Y-YY)  703.22.22 

22  YV«Y 
703  CONTINUE 

CALL  MATPRO(LOCM.A.¥C,PROO) 

00  370  IPM.LOEH 
00  370  JPM.LOEN 

Y*FL0A7(J-1).ZZ(IP.JP1‘RELLU1»VC(I».JP»-CNPXX(I)»VR(IP.JP) 

Y.Y-PROO(IP.JP) 

YaOABS(Y) 

23  IF(Y-YY)  370.24.24 

24  YY»Y 

370  CONTINUE 
717  00  307  IPM.LOEN 
00  307  JPM.LOEH 
VR(IP.  JPI-ZUP.  J>) 

¥C(IP.JP1«ZZ<IP.JP) 

307  CONTINUE 

IF(J.NE.MHM)  GO  70  149 
CALL  NATPRO (L0E4. A. VM.PROO) 

00  893  IPM.LXH 
00  843  JPM.lOEN 

Y»HELL(I)»VR(IP.JP)  *CHPXX(I>4¥C(IP.J»)*PROO(IP.JP) 

Y«OA»S(Y) 

41  IF(Y-YT)  *93.42.42 

42  YY«Y 
893  CONTINUE 

CALL  NATPROaOCN.A.VC.PROD) 

00  943  IPM.LOEN 
00  993  JP«1.L0EN 

Y«CNPXX(l)»¥R(I».JP)-RELL(I>»VC(lP.JP).PROOtIP.JP» 

Y«0A8S(Y) 

43  IF(Y-YY)  993.44.44 

44  YY«Y 
993  CONTINUE 
109  CONTINUE 
809  CONTINUE 

RRI TE (4.7557 
WITE(*.19)  YY 

19  FORNAT ( 3X . 30NNAK I NUM  ENTRY  OF  ERROR  MATRIX*. 1015.0) 

GO  TO  066 

179  FORNAT  (///3X.3NZ7 (.12. IN. .12*1  IN)  MATRIX  IS/) 

77*  FORNAT (///3X.2HZ( . I2.1M. .12.11 N)  MATRIX  IS/) 

1 F0RNATU2) 

76*  CONTINUE 
ENO 
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SUBROUTINE  5INVRT(A«W«0»L.M) ’ 

c 

C THIS  raUTINE  CALCULATES  the  INVERSE  or  * HATRTJC 

c 

5 OINENSION  A(I).L(I).M(I)  - 

DOUBLE  PRECISION  A.8IGA.0.HOL0  

NK»“N 

00  80  K*I.N  

10  NK=NK*N  ' > 

“ L(K)*K  ' ■ ' “ — ...  — : 

M(K)=S 

......  -kx»NK«K'  — ' 

BIGA*A(KK) 

15  00  20  J«K,N  : 

IZ*N«(J-I) 

DO  20  I*K.N  ' ~ - - 

ijmz*i 

10  IF  (DABS  (BIGA)-DABS  <*"1 IJ)  )T  15.20.20  

20  15  BI6A«A(IJ) 

LOOM  - - 

M(K)*J 

20  CONTINUE  ' ~ ' 

C INTERCHANGE  RONS 

>ST  • " • J*L(K)  • ■ • • — 

IF(J-K)  35.35,25 

25  KI*K-N  . . ..  ..... 

00  30  IM.M 

KI*KI»N  

TO  HOLO»-A(KI) 

JI.K1-K.J  - 

A(KI)*AMI) 

30  A(JI)*HOLD  

C INTERCHANGE  COLUMNS 

15  35  I«M(K»  - 

IF(I-K)  A5.A5.38 

38  JP*N«(I-tr  

DO  AO  JM.N 

JK*NK*J  

iO  JI»JP»J 

MOLD*-A<JX) 

A(JK)*A(J1) 

• — ait  ATjt»«HOLo  — - - - 

C DIVIDE  COLUMN  8»  MINUS  PIVOT  (VALUE  OF  PIVOT  ELEMENT  1$  CONTAINED 

AS  C IN  BIOAI 

A5  IF(BIGA)  A9.A6.A8 

A6  0*0.  - 

RETURN 

A8  DO  SS  IM.N  — 

SO  IF(I-K)  50.55.50 

SO  IK*NK«I 

A ( IK)  *A  ( IK)  /( >8 IGA) 

55  CONTINUE  - - ■ • 

C REDUCE  MATRIX 

55  DO  65  IM.M  - 

IK*NKM 

HOLO*A(IK) 


Figure  5.  (Continued). 
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OPTsl 


60 


65 


75 


90 


^5 


loo 


60 

62 


65 

C 


70 

75 

C 


00 

C 

100 

10S 

ioe 


110 

120 

125- 


MO 

ISO 


"fJPl-N  --  - - - - 

00  65  J*1,m 
“ lj=IJ*N 

IF(I-K)  60.65,60 
IF(J-K)  62*65*62 

kj=ij-i*k 

c on  n nue  LT>  * 4<  KJ  r*  * ' rj  

DIVIDE  row  3Y  pivot 
kj=k-n 
OD  75  J=1,N 
KJsKJ.N 

IF(J-K)  70,75.70 

a<kj>*a<kj) /bioa 
CONTINUE 

PPODUCTa  OF  PIVOTS 
0=D*dIGA 

er  ,ec!pro«l 

CONTINUE 

final  row  ano  column  interchange 

K--(K-i) 

IF(K)  150,150,105 
I=L  K) 

IF(I-K)  120,120.108 

UQ=N* (K-l ) 

JR=N»(I-l) 

00  110  J=1,N 
,"<=JO*J 

- Tt0L0=A(JK)  - _ __ 

JI=JR»J 
A ( JK)  *-A  ( JI ) 

A{JI)»HOLO 

J=M(K) 

IF(J-K)  100,100,125 
l\l=K-N 

00  130  1*1, N 

KI=KI*N 

ROLO*A(kI) 

JI=KI-K.j 

A(KI)*-AUI) 

*(Jl)«H8tO 

GO  TO  100 
RETURN 

End 
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tu 


SUBROUTINE  ARRAY  tMOD£-»T.-UiMfN.S.D)  ' 

C 

C THIS  ROUTINE  PREPARES  A MATRIX  FOR  SINVERT 
C 

oimension  siiftOdi 

DOUBLE  PRECISION  S.O 

— 

C TEST  TYPE  OF  CONVERSION 

ircMOor-dr  100.1001120  ~ ■ 

C CONVERT  FROM  SINGLE  TO  DOUBLE  DIMENSION 
100  IJbI»J*1 
NM-N*J*l 

00  I10  K=1VU 

NM*NM-NI 
DO  110  LS1 tl 
IJ=IJ-1 
NM-NM-1 

110  0 (NM>  =S ( I J) 

GO  T0“T40  

C CONVERT  FROM  DOUBLE  TO  SINGLE 
120  I J=0 
NM=0 

00  130  K»l,J 
00  125  L=l.I 

IJ=IJ*1  - ' 

NMsNM* 1 

125  S(IJ)*D<NMJ 
130  NM=NM*NI 
140.  RETURN 
END 
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SUBROUTINE  MATPR0(L0EM,A,8.PR0D) 

C 

C THIS  ROUTINE  CALCULATES  THE  PRODUCT  OF  TWO  MATRICES 

C 

OIMENSION  A tl 0. 10) *8(10*10) »PROO (20.20) 

OOUBLE  PRECISION  A. 8.  PROD 

DO  20  m.LDEM  ' 

DO  20  J=1.LDEM 
PROOIT.JMO. 

00  20  L«1.L0EM 

20  PROO(I.J)=PROD(I.J)*Atl.L)*B(L.J) 

RETURN 

END 
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SUBROUTINE  WRITE1  (IQ. A) 

C 

C THIS  ROUTINE  «RITES  A MATRIX 

c 

S DIMENSION  A ( 1 0 * 1 0 ) 

OOUBLE  PRECISION  A 

00  1 1*1.10 

1 WR1 TE (6.2)  (A(l.J) ,J=1,1Q) 

2 FORMAT (6I4X.D16.9) ) 

10  RETURN 

ENO 


Figure  5.  (Continued). 
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SUBROUTINE  COMCALfAVB.C.O.E.F)  - 

C 

c this  routine  calculates  the  product  of  a double  precision  complex 

C NUMBER 
C 

OOUBLE  PRECISION  A.8.C.0.E.F 

E=A«C»B»0-  — - 

F=A#0*8#C 

RETURN  

ENO 


SUBROUTINE  VANDR  74/74  OPT*l 


FTN  4.2*75067  05/16/75  10.02.53. 


SUBROUTINE  VANOR (RELL.N.M.K, VI  ‘ 

C 

C THIS  ROUTINE  CALCULATES  THE  VANOERMONDE  FOR  A MATRIX  KITH  ONLY  HEAL ' 

C EIGENVALUES 

5 C 

OIMENSION  RELL(I)*M(l)tV(20*20) 

OOUBLE  PRECISION  REEL, V " “ 

OOUBLE  PRECISION  Z.ZZ 

NSUM=0 

]0  00  10  L=1 «K 

MM=H(L1 

IF(L.EQ.l)  NC*1 
IF(L.GT.l)  NC-NSUM*T 
NSUM-NSUM*MM 

15  ZZ=1.00 

Z=RELL (L> 

00  30  1*2. N 

zz=zz*z 

Vd.NCI=ZZ 

?0  30  CONTINUE 

vii.Ncm.oo 
IF(MM.EQ.l)  50  TO  10 
00  40  1=1. N 

00  40  U=2.MM 

35  ■ JJ=NC*J-1 

IF(1-J1  97,98.99 

97  V(I,JU)=O.DO 
GO  TO  40 

98  V(1.UJ1=I.00 

10  GO  TO  40 

_ 99  v(i,uji=v(r-T.jj-n*vfi-i,ujr*z  _ ■ 

40  CONTINUE 
10  CONTINUE 
RETURN 

15  ENO 


SUBROUTINE  NRITEP  74/74  0PT=1 


FTN  4.2*75067  05/16/75  10.02.58. 


SUBROUTINE  VRITEP(XCOF.N) 

C 

C THIS  ROUTINE  WRITES  A POLTNOMlAL 

C 

5 UIMENS10N  XCOF(l) 

OOUBLE  PRECISION  XCOF 

..  -TOnert.nm-xcomi  ■■■■•■ 

NN«N-I 

DO  10  1=1, NN 

10  10  WRITE (6.1001  XCOFd.il, I 

WRITE (6*101)  XCOFIN.ll.N 

100  FORMAT </3X,lH(, 1015. 8.IH1. IX, 4H=X»».12,1X,IH*1 
- 10rT0RMArr/3X,IH(Vrai5.8,IH>,IX,4H«X»o,12) 

103  FORMAT (/3X,IH(, 1315. 8. 1H), IX, 1H*1 
15  RETURN 

ENO 


at 


Figure  5.  (Continued), 
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SUBROUTINE  SORT 


74/74  OPT=l 


FTN  4.2*75067 


05/16/75  10.03.01. 


10 


15 


70 


">0 


15 


iO 


SUBROUT  TNE-SORT  (ROOTR. ROOTI.  RELL,  CM°XX,  STORX, ST  ORT'.M-.N.IC.OIO 


THIS  ROUTINE-SORTS  THE  ROOTS 


40 


DOUBLE  PRECISION  OIC 
OIMENSION  Mil) 

TIIMENSroN” STTOXTT)  .STORY (I ) — 

OIMENSION  ROOTR ( 1 ) .ROOT I ( 1 ) . RELL <} > .CMPXX ( 1 ) 
0CU8LE  PRECISION  ROOTR, ROOT! . RELL. CNPXX 
OOUBLE  PRECISION  STORX. STORY 
00U8LE  PRECISION  X, Y.rf.Z.00,06 
K=0 

KK=0  — • 

DO  10  1*1, N 
KrK*  1 
M(K)=0 
X*ROOTR ( 1) 
r*ROOTI(l) 

RELL (K) =X  ' '•  

CMPXX <K)=Y 

NN=N-KK 

NK=0 

00  20  J=1.NN 

w=rootr( j) 

Z=ROOTI(J) 

oo=oabs(x-ri 

08=DA8S(Y-Z) 

ir(DO.LT.OK.ANO.ja.LT.OK)  GO  TO  30 

NK=NK*1 

STORX (NK)*W 

STORY <NK)=Z  '■ 

GO  TO  20 

M(K)=M(K)*1 

CONTINUE 

00  40  L*1,NK 

R00TR(L)*ST0RX(L) 

ROOTI (L) *ST0RY  fL) 

CONTINUE 


KK*KK*M(K) 
IF(KK.EO.N) 
10  CONTINUE 
SO  RETURN 
END-  - 


GO  TO  50 


Figure  5,  (Continued), 
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SUBROUTINE  POLRT 


74/74  0PT=1 


FTN  4*2*75067 


05/16/75  10.03.04. 


10 


15 


?0 


25 


30 


35 


.0 


45 


50 


55 


SUBROUTINE  -POtRTtXCOF  .COf  .N,  BOOTS, crOOT  t I£R) 

this  routine  calculates  the  roots  of  a polynomial 

OOUBtE  PRECISION  FT  - - 

OOUBLE  PRECISION  XCOF.COF.ROOTR.ROOTl 
-OOUBLE  PRECISION  XO.TO.r.-T.-XPR.TPR 
OOUBLE  PRECISION  OX. OY, TEMP, ALPHA 
DIMENSION  XCOrrtr.COF  U ) VROOTR  a ) .ROOTl  ( 1) 

1F1T*0 

N*M ' 

1ER*0 

— TRXC0FCN*irnffV7SVIff  - 
10  IF(N)  15.15,32 
15  1ER=1  - ----- 
20  RETURN 

25  I£R*4  ' - - 

GO  TO  20 

30  IER=2 

GO  TO  20 

3 Z IF1N-361  35.3S.30 
35  NX=N 

NXX=N*I' 

N?=l 

- -KjT*7r*i 

00  40  L=1,KJI 
MT=KJI-L»I 
40  COF(MT)*XCOF(U 
45  XO=. 00500101 
YO*. 01000101 

50  X=XO 

XO=-IO.O»YO 
YO="l0.O*X 
X=XO  — - 
Y*YO 

INslN»T - — ' 

GO  TO  59 
55  IF1T=C 
XPR=X 
YPR*Y 

59  1CT=0 

60  UXSffflO 

UY=0.0 
V-0.0 
YT=0.0 
XT-I .0 
U=COF<N*I) 

IF  CUT  6S,  130703 

65  DO  70  1*1, N 
L*N-1*I 
TEMPsCOF(L) 

XT2*X»XT-Y»YT 
YT2»X»YT*Y»XT 

U*0*-TEMP*XT2-  

Y*Y*TEMP»YT2 
FI»I 


j 


Figure  5.  (Continued). 
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SUBROUTINE  POLRT 


fin  <*.2*7SQ6V 


OS/16/75  ^/U. 


UX*UX*FI*XT»TEMP 
uy*uy-fi*yt»temp 
XT*XT2 
70  YT*YT2 

sumsg*ux*ux*uy*ut 

IF  (SUHSQ)  75*110*75 
75  DX*(V*UY-U«UX)  /SUHSQ 
X*X*OX 

1>Y*-(U*UY*V*UX) /SUHSQ 
Y*Y*OY  v 

70  IF(0ABS(0Y)*DABS(0X)-1.0-12)  100.B0.B0 
80  ICT*1CT*1 

IF ( ICT~5QG ) 60.BS.BS 
85  1F<IFIT)  100*90*100 
90  IF ( IN-5)  50*95. 9S 
95  l£R*3 
GO  TO  20 

100  00  10S  1*1 .NXX 

mt=kji-l*i 

T£MP*XCOF(MT) 

XCOF(MT)*COF(L) 

105  COF (L) *T£MP 
ITEHP*n 
N*NX 

NX* I TEMP 

IFUFIT)  120.S5.120 
110  IF(!F1T)  1 IS.SO * 1 IS 
1 IS  X*XPR 

y*ypr 

120  1FIT=0 

122  IF (OABS ( Y) -1 .0-1 0*0 ABS ( X) ) 13S.125.12S 
125  ALPHA*X*X 

sumsq*x*x*y*y 

N*N-2 
GO  TO  140 
130  X*0 • 0 
NX*NX-1 
NXX*NXX-1 
135  Y*0.0 

SUHSQ* 0.0 

ALPMA-X 

N*N-1 

1 AO  C0F(2)*C0F(2)*AlPHA*C0F(1) 

145  00  150  L*2.N 

150  COF  (L* 1 ) *COF (L* 1 ) *ALPMA*COF (L) -SUMS3*C0F (L-l ) 
1SS  R00TI(N2)*Y 
ROOT* (N2)*X 
N2*N2*1 

IF (SUHSQ)  160.165.160 
160  Y»-Y 

SUMSU*0.0 
GO  TO  1S5 
165  IF (N)  20.20* 4S 
ENO 


SUBROUTINE  TRACE 


74/74  0PT*1 


FTN  4.2*7S067 


05/16/75  10.10.24. 


Subroutine  trace(a.n.tr) 

THIS  ROUTINE  CALCULATES  THE  TRACE  OF  A MA1R1X 

01 MENS ION  4(10.10) 

OOUBLE  PRECISION  SUM.A.TR 

sum *0.00 
00  10  1*1. N 
10  SUH*SUH*A (1.1) 

TR*SUH 

RETURN 

ENO 


SUBROUTINE  TRACE2  74/74  0PT*1 


SUBROUTINE  TRACE2 < A.N. TR) 


FTN  4.2*7S067 


OS/16/75  10.10.29. 


THIS  ROUTINE  CALCULATES  THE  TRACE  OF  A MATRIX 

DIMENSION  A (20.20) 

OOUBLE  PRECISION  SUM.A.TR 
SUM*0«l)0 
00  10  1*1 »N 
10  SUM*SUH*A(I.I) 

TRsSUM 

RETURN 

ENO 


Figure  5.  (Concluded), 
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